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FREE-FORM AIRFOIL SHAPE OPTIMIZATION UNDER UNCERTAINTY USING 
MAXIMUM EXPECTED VALUE AND SECOND-ORDER SECOND-MOMENT 

STRATEGIES 


LUC HUYSE* 

Abstract. Free-form shape optimization of airfoils poses unexpected difficulties. Practical experience 
has indicated that a deterministic optimization for discrete operating conditions can result in dramatically 
inferior performance when the actual operating conditions are different from the - somewhat arbitrary - 
design values used for the optimization. Extensions to multi-point optimization have proven unable to 
adequately remedy this problem of “localized optimization” near the sampled operating conditions. This 
paper presents an intrinsically statistical approach and demonstrates how the shortcomings of multi-point 
optimization with respect to “localized optimization” can be overcome. The practical examples also reveal 
how the relative likelihood of each of the operating conditions is automatically taken into consideration 
during the optimization process. This is a key advantage over the use of multipoint methods. 

Key words, airfoil shape optimization, sensitivity analysis, statistical decision making, robust design, 
stochastic optimization, second moment analysis, expected value optimization, multi-point optimization 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. The design process of a structure or device is essentially a decision-making process. 
Appropriate values of design variables, which optimize the performance of the design, need to be selected. The 
specification of one or more design operating conditions allows the engineer to use deterministic optimization 
schemes. An example hereof is found in airfoil design where a cruise Mach number and target lift coefficient 
are specified and the objective is to minimize the drag under those constraints. 

Optimization of an analytical model is a tool to develop better designs. Recent advances in computing 
power and the development of more accurate computational fluid dynamics codes (CFD) should, at least in 
theory, allow one to compute the optimal shape of an airfoil for a particular application. Unfortunately, the 
use of deterministic optimization techniques leads to unexpected problems and often unacceptable results. 

An important concern in the shape optimization of airfoils is the sensitivity of the final optimal design 
to small manufacturing errors or fluctuations in the operating conditions. Tightening the tolerances in 
the manufacturing process may prove prohibitively expensive or practically impossible to achieve. It is 
expensive to produce such a precision design and impossible to maintain this pristine shape during routine 
flight operations. Moreover, a certain variability in the operating conditions (e.g. cruise Mach number) 
cannot be avoided. The sensitivity of the design performance to relatively small uncertainties provides an 
incentive to use so-called “robust” optimization methods, which directly assess the effects of the uncertainties 
on the design performance. 

Several different non-deterministic approaches (Taguchi methods, bounds-based, minimax, fuzzy and 
probabilistic methods) can be taken to achieve “robustness”. This research deals with robust design of airfoils. 
Other researchers have worked in the field of reliability-based design, and have demonstrated optimization of 
aerospace structures (e.g. [7]). Figure 1.1 explains the difference between robust and reliability-based design 
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research was supported by the National Aeronautics and Space Administration under NASA Contract No. NAS1-97046 while 
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Fig. 1.1. Uncertainty classification 


problems. In an engineering context, risk is typically defined as the product of the probability and the impact 
or consequence of an event [16]. Structural reliability techniques are particularly useful to assess the risk 
associated with infrequent but potentially catastrophic events (such as turbine blade failures [1]). On the 
other hand, robust optimization techniques account for the impact of everyday fluctuations of parameters 
(such as variations in operating speed) on the overall design performance assuming that no catastrophic 
failures occur. 

2. Achieving “Robustness”. Developing optimization methods which result in more “robust” designs 
sounds appealing. The term “robustness” has been coined to mean a variety of things. This section outlines 
some of the more popular goals of “robust” design. Loosely speaking we can distinguish the following 
meanings: 

1. Identify designs which minimize the variability of the performance under uncertain operating con- 
ditions. This is the objective of Taguchi methods [12]; they are most practical when the mean value 
of the performance can be adjusted at negligible cost. 

2. Mitigate the detrimental effects of the worst-case performance. This is the objective of MiniMax 
strategies [23]. They choose the design with the “best” worst-case performance. 

3. Provide the best overall performance over the entire lifetime of the structure or device [20]. This is 
the objective of this work. 

In the remainder of this paper the focus is on the third definition: robust optimization which tries to 
achieve the best performance (or minimal cost) for all possible combinations of the operating conditions. The 
paper focuses on the effectiveness of the optimization strategy rather than on particular implementations of 
optimization algorithms. 
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To our knowledge, non-deterministic approaches are quite new to aerodynamic optimization. First, 
an overview of existing deterministic attempts at introducing robustness is presented. Subsequently, we 
introduce an inherently statistical approach based on Van Morgenstern’s Maximum Expected Value Criterion 
[20]. Numerical results for a 2D airfoil shape optimization problem obtained using the different optimization 
strategies are compared. 

3. Airfoil Optimization Problem in Transonic Regime. 

3.1. Problem Formulation. In this section we present a transonic airfoil optimization problem to 
which the various optimization strategies will be applied. All optimizations use the NACA-0012 profile 
as the baseline geometry. The NACA-0012 is splined using 23 control nodes. The design variables in the 
optimization problem are given by the vertical positions of the control nodes and the angle of attack a. Three 
control nodes are in locked positions: one at the leading edge, and a double control node at the trailing edge. 
Box constraints limit the maximum movement of each of the 20 spline control nodes (see Figure 3.1). The 
inviscid Euler equations for the flow are discretized on unstructured meshes [4]. The sensitivities of both lift 
and drag with respect to the design parameters are efficiently calculated using a discrete adjoint formulation 
[3]. 



Fig. 3.1. Design variables and box constraints for the airfoil 

The objective is lift-constrained ( Cy > Cf) minimization of the drag Cd over the Mach range M e 
[0.7, 0.8]. 


rnindeD C<y(d,M) over M € [0.7, 0.8] 
subject to Ci > 


where d is the vector of design variables and D is the design space. In this study, the Mach number is 
the only uncertain operating variable; no additional model uncertainties are included. No constraints are 
imposed on the pitching moment C m . Because the NACA-0012 is not suitable for the transonic regime, 
substantial improvements are to be expected, especially when the target lift Cf is relatively high. 

3.2. Grid Convergence. Four different grid resolutions were used to check the convergence of the 
flow solution, especially the drag (7<y. The far field boundary for grid 1 (Figure 3.2) is located at 10 chord 
lengths. For grids 2, 3 and 4 (see Figures 3.3, 3.4, and 3.5) the far field boundary is located at 50 chord 
lengths. For each of the grids the number of elements along both the airfoil inner boundary and the circular 
far field boundary are specified in Table 1. Grid generator details are given in [18]. Figure 3.6 shows the 
drag-profile in the Mach-range M e [0.7, 0.8] for a constant angle of attack a = 1.5°. Grids land 2 are 
considered low-fidelity models, while grid 4 is the high-fidelity or “truth” model. To limit the computational 
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cost, only grids 1 and 2 will be used for the actual optimizations. Upon completion of the optimization a 
check using grids 3 or 4 can be performed. The idea is to use the fast low-fidelity grids as much as possible, 
but to verify the results using the high-fidelity grids [2]. 

Table 3.1 

Grids used in examples (Grid 1 is generated using Double9, all other grids are generated using AFLR [18]) 



Airfoil elem. 

Far field elem. 

Nodes 

Faces 

Elements 

Grid 1 

129 

32 

411 

1170 

822 

Grid 2 

124 

32 

3060 

9025 

6120 

Grid 3 

250 

64 

9246 

27424 

18492 

Grid 4 

500 

128 

32067 

95574 

64134 


O 

> 



Fig. 3.2. Baseline NACA-0012 Airfoil: Grid 1 


Figure 3.6 reveals that grids 1 and 2 overestimate the drag. The overestimation of grid 1 has a double 
cause: it is an extremely coarse grid and the far-field boundary is located at 10 chord lengths only. The 
overestimation in grid 2 is less than for grid 1 and occurs in equal amounts irrespective of the Mach number. 
There is virtually no difference in the drag profiles computed on grids 3 and 4. Very similar grid-induced 
differences in the drag profiles were observed in the Mach sweep curves ( Cd as a function of M) for the 
optimized airfoils. 

3.3. Constraint Handling. In this problem we are not concerned with rapid Mach number variations. 
Only slowly varying Mach numbers (steady states) are considered. Because the Mach number is stable, the 
angle of attack can be adjusted to reach the required lift Cf. Consequently, the constraint in Eq.(3.1) is not 
probabilistic but remains deterministic. If point-in-time fluctuations are considered, the angle of attack can 
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Fig. 3.4. Baseline NACA-0012 Airfoil: Grid 3 
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no longer be adjusted to assure that enough lift is provided. The constraint needs to be transformed into a 
probabilistic constraint, which can only be met with a specified probability. 

A mathematically straightforward approach to the optimization problem is to account for the inequality 
constraint in Eq.(3.1) using a feasible direction search algorithm. For this purpose, we linked FUN2D [4] to 
several optimizers - Modified Method of Feasible Directions (CONMIN) and Sequential Linear Programming 
(ADS/SLP) - in iSIGHT v.5.5, a product of Engineous Software Inc. [11]. However, the lift constraint can 
also be eliminated if we take advantage of the specific nature of the problem. We solved the resulting 
unconstrained optimization problem using PORT, a bound constrained trust region algorithm [14]. 

4. Deterministic Approach to Airfoil Shape Optimization. 

4.1. Formulation of Deterministic Optimization Problem: Single-Point Optimization. In 

a deterministic context, aerodynamic shape optimization of airfoils is concerned with obtaining the most 
aerodynamically favorable geometry for fixed - either known or assumed - operating or design conditions. 
Consider the practical case where the drag Cd is to be minimized at a given, fixed free flow Mach number 
M: 


min deD C d (d,M) 

subject to Ci,i( d, M) > Cf for i = 1, ..., n 


This deterministic, single-point optimization model is not necessarily an accurate reflection of the re- 
ality. The formulation in Eq.(4.1) contains no information regarding off-design condition performance. It 
is documented by other researchers [9] that, with formulation Eq.(4.1), the drag reduction is attained only 
over a narrow range of Mach numbers (see Figure 4.1a.). We will refer to this in the remainder as ’’localized 
optimization”. Drela explains that the optimizer creates a “bump” on the airfoil to fill the transitional 
separation bubble (see Figure 4.1b). This effectively reduces the drag penalty which occurs when a bubble 
undergoes transition and reattachment [8]. However, the location of this bubble varies with M and this 
explains the poor behavior in off-design conditions for this “locally optimized” design. 




Fig. 4.1. Single point optimization using Grid 2 and Cf = 0.2: drag profile and airfoil geometry 


It can be concluded that the real problem is not with the optimization code, which is likely to perform 
just fine, but with the problem formulation of Eq.(4.1). The “local” optimization effect is particularly 
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worrisome if substantial variability is associated with the operating conditions. Trade-offs between different 
design conditions should explicitly be considered in the problem formulation. 

4.2. Multi-Point Optimization. A straightforward, but heuristic, approach to avoid point-optimized 
designs is to consider different Mach numbers and to generalize the objective in Eq.(4.1) to a linear combi- 
nation of flight conditions (m in total): 


min deD WiC d {d, Mi) 

subject to Cij( d, M , ) > Cj) for j = 1, ..., n 


Practical problems arise with the selection of the flight conditions M* and with the specification of the 
weights Wi . There is no clear theoretical principles to guide the selection, which is in fact largely left up to 
the designer’s discretion (see e.g. Refs. 6, 9 and 10). 

With the multi-point formulation of Eq.(4.2), an improved Cj can be realized over a wider range of 
Mach numbers M [9]. However, this formulation is still unable to provide a truly global solution by avoiding 
localized optimization. In fact multiple “bumps” appear on the airfoil, one associated with each flight 
condition Mj. In the transonic regime, each bump occurs at the shock foot location for each of the sampled 
Mach numbers. 

A clear example of “localized optimization” is given in Figure 4.2. This optimization was performed on 
Grid 2, with lift constraint C* = 0.4, using equal weight (w* = 0.25) for each design condition. Because the 
drag is highest at the higher Mach numbers, the drag at those Mach numbers has been reduced significantly 
during the optimization process. Figure 4.2 clearly indicates that these drag improvements occur only at 
particular Mach numbers and are rapidly lost when the actual Mach number deviates from any of the selected 
design Mach numbers. 

Figure 4.3 explains this in more detail using 2 contour plots of the local Mach number. The operating 
conditions (free flow Mach numbers) are very similar but the flow solutions (local Mach number) are very 
different. The multi-point optimization process introduces geometric features to the airfoil which lock the 
shock waves in place. Since we used four point optimization, we have 4 shocks (see also [9]). In the figure 
on the left, four shocks can be distinguished along the top surface. In the figure on the right, the most-aft 
shock has basically disappeared. The optimizer has locally modified the geometry and eliminated the shock 
associated with Mach 0.8, which was one of the design conditions. 

4.3. Heuristic Approach: Weighted Average of Geometries. Since derivatives have to be calcu- 
lated for each variable at various operating conditions, the computational cost associated with multi-point 
optimization using CFD is substantial. Knowledge-based methods, which incorporate some “common-sense 
engineering knowledge” into the optimization process, may be a valuable and cheaper alternative. One such 
example is the method of the Weighted Average of Geometries (WAG) [6]. 

In this method the optimal design is obtained as a weighted average of n single-point optimal designs, 
each one of them corresponding to one of the n chosen operating conditions. The weighting factors depend 
on the relative importance of each operating condition (similar to the wfs in Eq.(4.2) in multipoint design) 
and on the assumptions made regarding the variation of the drag with geometry and Mach number (response 
surface model). 

The method is computationally efficient, but suffers some of the same drawbacks as the multi-point 
optimization: 
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Fig. 4.2. Drag profile for a four-point optimization using Grid 2 and Cf = 0.4, sample points are M = 0.7,0.733,0.766 
and M = 0.8 


Operating condition: free flow Mach number is 0.795 Operating condition: free flow Mach number is 0.8 
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Fig. 4.3. Local Mach number for two free flow conditions after 4-point optimization using Grid 2 and C ; * = 0.4 


1. As with multi-point optimization, the selection of design conditions to be included in the optimization 
process is arbitrary. 

2. Single point optimizations result in bumpy airfoils. Linear combinations of bumpy airfoils will remain 
bumpy. 

3. Campbell [6] assumes a quadratic variation of Ca with the geometry. This assumption may become 
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questionable when the geometries are much different due to large variations in Mach number. 

5. Noil-Deterministic Approaches. 

5.1. Design as a Decision-Making Process. Designing a structure or device is essentially a decision- 
making process. Appropriate values of the design variables d need to be selected which optimize the perfor- 
mance or the utility of a design. The designer has full control over the design variables, such as the geometry 
of the structure and the type and grade of materials used for it, but the operating conditions of a structure 
or device, such as the loads or the operating speeds, will typically vary during the design life time. 

Since each operating condition parameter may take on a range of values over the lifetime of the design, it 
is possible to collect their histograms (and joint histograms). From a subjectivist point of view, the operating 
conditions 8 are then effectively modeled as random variables. 

The previous section indicated that a specific design may perform exceptionally well for a selected 
operating condition, say the free flow Mach number M, but may perform poorly for slightly different values 
of M, which are quite likely to occur. The impact of the uncertainty of M on the design performance should 
be taken into account when the quality of a particular design is assessed. 

5.2. Bounds-Based Methods. Bounds-based methods recognize that a designer does not have access 
to data with unlimited precision: some of these data are intrinsically variable (operating conditions), some 
of them can be collected or measured with limited precision only (sampling error), and sometimes the 
measurement process itself is imperfect which introduces bias. Gu et al. [13] recognize that analytical and 
numerical models fail to yield the “correct” result; they call the difference between the average experimental 
data and the value obtained from the analysis code the bias error. All other errors are lumped into “random 
variability” . They assign bounds to each of those uncertainties and develop a methodology to incorporate 
the bounds directly into the multi-disciplinary optimization [5]. 

In our opinion, the main drawback of bounds-based methods is that they essentially assume that all 
outcomes within the specified bounds are equally likely to occur. This conflicts with the intuitive sense that, 
by their very definition, extremes should occur with much lower frequency than average or “normal” behavior. 
This implies that error bounds will grow quite rapidly as more and more uncertainties are explicitly considered 
in the design or optimization process. An explicitly statistical approach takes the relative likelihood of 
combined extremes versus other joint occurrences into account by means of the Probability Density Function 
(PDF). In short two problems can be associated with the bounds-based approach: 

1. Bounds cannot always be identified accurately. 

2. Bounds-based methods are over-conservative. 

5.3. Explicitly Statistical Problem Formulation. Reconsider the basic problem in Eq.(3.1): min- 
imize the drag Cd over a range of free flow Mach numbers M while maintaining lift Ci > Cf . Note that M 
is now treated as a random variable. The optimization problem Eq.(4.1) is now interpreted as a statistical 
decision-making problem. 

In the presence of uncertainty a designer is forced, in effect, to take a gamble. Under such circumstances, 
rather than naively hoping for the best or over-conservatively focusing on the worst, the right decision consists 
in the best possible choice of the design, whether favorable or unfavorable operating conditions occur. All 
decision problems have two essential characteristics [20]: 

1. A choice, or sequence of choices, must be made among various possible designs. 

2. Each of these choices corresponds to a performance, but the designer cannot be sure a priori what 
this performance will be. The exact performance also depends in part on unpredictable events, in 
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Fig. 5.1. Decision tree example 


this case the operating conditions. 

In our example, only one initial choice (i.e. airfoil selection) regarding the design needs to be made (see 
Figure 5.1). However, the method is easily extended to accommodate sequential decisions. For instance, we 
may want to run some wind tunnel experiments in the hope of reducing the uncertainty on some mathematical 
model predictions. A decision tree would indicate if the upfront cost of the prototypes is offset by the expected 
reduction in parameter uncertainty. 

According to the Von Neumann-Morgenstern statistical decision theory [20], the best course of action in 
the presence of uncertainty is to select the design which leads to the lowest expected drag. This is commonly 
known as the Maximum Expected Value criterion (MEV). The risk p , associated with a particular design 
d, is identified as the expected value of the perceived loss associated with the design. The best design or 
decision, which minimizes the overall risk, is referred to as the “Bayes’ decision” . In our problem formulation, 
the Bayes’ risk p* and Bayes’ decision d* are given bu (5.1) and (5.2) respectively: 


p* = min d(E £> J M Cd(d, M)f M (M)dM 
subject to Ci (d, M) > Cf 

P*= f C d {d*,M)f M (M)dM 
J M 

where / m(M ) is the probability density function of the free flow Mach number M. 

The practical problem with formulation Eq.(5.1) is that integration is required in each of the optimization 
steps. Since the objective function C d is computationally expensive to evaluate, this approach, although 
theoretically sound, becomes prohibitively expensive. Therefore a computational scheme that minimizes the 
number of function calls is desirable. 


(5.1) 
c 

(5.2) 


li 



5.4. Analytic Approximation of the Expectation Integral. When the variability of the free flow 
Mach number M is not too large, a second-order Taylor series expansion of C d around the mean value M 
may be a sufficiently accurate model of the variation of the drag Cd with respect to M. 


(5.3) Cd{ d, M ) = C d ( d, M ) + Va f C d .(M - M) + | V 2 M C d .(M - M) 2 

When substituted in the Bayes’ risk expression (5.1), the linear term VMC d .(M — M) in Eq.(5.3) 
disappears after integration over M because the Taylor series is built around the mean value M. The Bayes’ 
risk is: 


(5.4) 


p = nun de £> 


C d (d,M) + Var(M) 


d,M 


subject to C'i (d, M ) > 6)* 


where Var(M) denotes the variance of the Mach number M. 

It seems that we have substituted an integration with an almost equally expensive computation of 
a second-order derivative. Even though the approximation may result in only moderate computational 
savings, this theoretical result provides additional insight into the problem. It follows from Eq.(5.4) that 
the variability of M can affect the optimal design only if the objective function C d is highly non-linear in 
this parameter. This is the case near the drag divergence Mach number Mdiv-i where the drag undergoes a 
sharp increase (recall Figure 4.1a). 

In mathematical terms, the advantage of working with expected utilities is that the minimum is second- 
order accurate with respect to variations in the parameters. This ensures a more global solution and “localized 
optimization” will be avoided. This can also be explained in an intuitive manner: the second-order derivative 
is a measure for the curvature. Since this curvature is now a part of the objective function, a design which 
results in a drag trough or “cusp” as found in the optimal solution in Figure 4.1 will not be accepted by the 
optimizer. The high curvature of the “cusp” would increase the objective in Eq.(5.4) and excessive localized 
optimization will be avoided. 

5.5. Direct Numerical Evaluation of Expectation and Comparison with Multi-Point Opti- 
mization. The integration with respect to M in Eq.(5.1) can also be performed numerically. Irrespective 
of the chosen integration scheme, the integral (5.1) can formally be written as (n k integration points): 


( n k 

w k .C d { d, M k ) + e(rifc) 

where the integration error e(n k ) — ¥ 0 as n k — ¥ oo. 

Formulation Eq.(5.5) is strikingly similar to Eq.(4.2). It is therefore interesting to analyze how the 
Bayes’ decision d* compares with the multi-point solution and exactly how localized optimization is avoided. 

In the multi-point approach the Mach numbers and weights need to be selected by the designer. In 
the statistical approach, the weights are directly related to the relative importance of each Mach number 
through the integration over the probability density. Which Mach numbers are used in the optimization 
depends on the chosen integration scheme. In short, the statistical approach removes the arbitrariness from 
the weighting process. 
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Careful comparison of Eq.(4.2) with Eq.(5.5) reveals the shortcoming in the multi-point formulation 
which causes localized optimization. Numerical integration of Eq.(5.1) results in Eq.(5.5) and includes a 
random, zero-mean error term e(n*), which decreases as the number of sampling points increases. The 
multi-point optimization Eq.(4.2) differs from Eq.(5.5) only in the sense that this error term is not explicitly 
considered in the objective function. However, omitting this error term from the optimization problem alters 
the structure of the problem at hand. The multi-point optimization effectively looks for the design, which 
minimizes the weighted sum of the goal function Cd , evaluated in the n* specified points M*. There is 
absolutely no control over what happens to the objective function Cd in the neighborhood around these n k 
sampling points. During the optimization iterations, the shape of the goal function Cd(d, M) is altered. As 
a result, the discrete sum in Eq.(4.2) may fail to be a good approximation of the integral in Eq.(5.1). 

In effect, multi-point optimization will prefer a design d, over a design d 2 even when design di is 
considerably worse than design d 2 in all but the specified sampling points. Multi-point optimization 
allows the optimizer to mold the goal function Cd to its own advantage. What was originally a random 
integration error is no longer random, and the discrete sum in Eq.(4.2) no longer approximates the integral 
in Eq.(5.1) at all. 

This rather annoying behavior is avoided if we can prevent the optimizer from exploiting the approxi- 
mation error in Eq.(5.5) to its own advantage. We need to make sure that the discrete sum in Eq.(5.5) really 
is an approximation of the integral in Eq.(5.1) at all times. In general terms, we need to ensure that the 
discrete sum Eq.(4.2) remains a good approximation of the integral in Eq.(5.1). An elegant solution is to 
randomize the sampling points M* in the evaluation of the integral but any adaptive optimization scheme 
that varies the location of the integration points Mk for each optimization step will do. Randomization of 
the integration points ensures that the optimizer maximizes the performance not just for n* specific values 
of A/*, but for any set of values M*, k =• -1, . . . n*. To minimize the loss of accuracy in the integration due to 
random location of the integration points, stratified sampling can be used to generate the Mk values. Our 
experience with the spline-based integration also suggests that the sampling points should not be allowed to 
be arbitrarily close to each other. 

5.6. Additional Considerations. The use of Eq.(5.5) for the optimization instead of Eq.(4.2) leads to 
numerical complications. Because of the random location of the integration points M*, a repeated evaluation 
of the objective function Cd for identical values of the design parameters d will lead to different results. This 
makes it hard to identify whether a new design is really better than a previous one, or if the “improvement” 
should be attributed to random fluctuations instead. When a trial solution d is still far away from the optimal 
solution d*, large improvements A Cd can be expected. This means that a very crude integration, which 
requires very few function evaluations, will suffice in the early stages of the optimization. The improvement 
of the goal function is expected to be smaller closer to the optimal solution, and more sampling points Mk 
will be required to keep the integration error small enough. Current research focuses on the development of 
a strategy which takes maximum advantage of this effect. 

In addition, the physical and mathematical models themselves used for the objective function will not 
be error-free. Each of these model errors can be treated as a random variable. Their effect on the optimal 
solution is readily assessed by extending the integration over these additional random variables. It is believed 
that the approximate second-order result in Eq. (5.4) will prove particularly useful for this purpose. In a first 
step we minimize the drag while keeping the model parameters fixed at their average level. If the second 
order derivative of the drag with respect to this model parameter at the solution of this simplified problem 
(without explicit error modeling) is sufficiently small, it can be concluded that these model uncertainties 
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will only have a minimal impact on the solution. Otherwise, the model uncretainties should explicitly be 
included in the problem formulation and a full integration is required. 

6. Application to a 2D- Airfoil in Transonic Regime. 

6.1. Various Optimization Strategies. In this section we compare results for the various optimiza- 
tion strategies presented in this paper. We apply the optimization techniques to the lift-constrained airfoil 
optimization problem introduced earlier. In this example we assume that the Mach number variations are 
bounded between 0.7 and 0.8. In particular, the following optimization formulations are compared: 

1. Optimization at a single Mach number: various Mach numbers are considered as the design point. 

2. Multipoint optimization using m Mach numbers. Each of the design conditions has pre-determined 
fixed weights Wi = 1/m, see Eq.(4.2). 

3. “Robust” optimization using a second-order second-moment approximation around the mean value 
M = 0.75, see Eq.(5.4). 

4. “Robust” optimization using stochastic integration. To allow for a direct comparison with the multi- 
point results, / m(M ) in Eq.(5.1) is a uniform distribution. Other distributions are considered as 
well, which demonstrates the versatility of the expected value optimization method. 

It is important to note that formulation 2 and 4 require a comparable computational effort. 

6.2. Single-Point Optimization Results. The single point case has 21 design variables: the angle 
of attack a and the vertical positions of the 10 spline control nodes at both the top and bottom surface of 
the airfoil. Figure 4.1a indicates a dramatic reduction of the drag Cj is obtained at M = 0.76 (also see the 
first entry in Table 6.1), but also reveals that this gain is rapidly lost when the free flow Mach number is 
away from this design value. 

Table 6.1 

Comparison of the drag reductions with respect to the original NACA-0012 at the design point and over the entire Mach 
range for single point optimization results. 


Optimization 

model 

Grid 

Target 
lift c ,* 

Design Mach 
number M 

A Cj, at the 
design point 

Actual total 
reduction A Cj 

CONMIN 

Grid 2 

0.175 

0.76 

44% 

12% 

PORT 

Grid 1 

0.175 

0.75 

24% 

11% 

SLP/ADS 

Grid 2 

0.2 

0.75 

49% 

21% 

SLP/ADS 

Grid 2 

0.2 

0.76 

60% 

18% 

SLP/ADS 

Grid 2 

0.4 

0.75 

84% 

52% 

SLP/ADS 

Grid 2 

0.6 

0.72 

85% 

30% 

SLP/ADS 

Grid 2 

0.6 

0.75 

87% 

53% 

SLP/ADS 

Grid 2 

0.6 

0.78 

91% 

78% 

SLP/ADS 

Grid 2 

0.6 

0.80 

75% 

57% 


The geometry plot in Figure 4. lb illustrates what happens. During the optimization a distinct “bump” 
is formed on the airfoil surface. The optimizer takes advantage of all degrees of freedom to achieve the 
lowest possible drag at the design Mach number M = 0.76, irrespective of what happens to the drag at 
other Mach numbers. Obviously, there is a penalty to be paid for this: even though the drag reduction at 
the design Mach number is 44%, the total average reduction over the entire Mach range is only 12% (results 
obtained with Grid 2 and ( 7 / = 0.175 using the CONMIN optimizer implemented in the iSIGHT package 
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At Mach 0.72 At Mach 0.78 



Fig. 6.1. Bump locations for different design Mach numbers for single point optimization (results obtained using Grid 2, 
C[ = 0.6, and SLP/ADS [22] optimizer) 



Fig. 6.2. Single point optimization results for various design Mach numbers (results obtained using Grid 2 , C [ = 0.6, 
and SLP/ADS [ 22 ] optimizer) 
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Fig. 6.3. Optimal drag profiles obtained using different 2-point optimization schemes (results obtained using Grid 2, 
wi = W 2 = 0.5, C ; * = 0.6 and SLP/ADS optimizer) 


[11]). The use of a single-point optimization would lead to the false conclusion that a 44% improvement has 
been achieved; the actual realized gain is only 12%. 

When unconstrained optimization using PORT [14] after elimination of the lift constraint is used (results 
for Grid 1 and C( = 0.175), the drag reduction at the design Mach number M = 0.75 is 24% and the actual 
reduction over the entire Mach range is found to be 11%. It can be concluded that both methods to handle 
the constraints lead to the same qualitative result. The numerical difference between the two results can be 
attributed to the overestimation of the drag of the baseline NACA-0012 in Grid 1, which results in a larger 
reduction in the optimized profile. 

Qualitatively similar results are obtained for other lift values and design Mach numbers and are summa- 
rized in Table 6.1. The localized optimization, illustrated in Figure 4.1, was previously documented by Drela 
[9]. Figure 6.1 shows that the bump on the airfoil moves aft when the Mach number increases. This pushes 
the shock towards the trailing edge and postpones the drag rise. Figure 6.2 indicates that the choice of the 
design Mach number considered in the optimization has an enormous impact on the final drag profile and 
airfoil geometry. The figure shows that simply selecting the average Mach number (so called Mean- Value 
Optimization at M = 0.75) or the Mach number with highest drag (M = 0.80) does not guarantee good 
overall performance over the entire Mach range. 

6.3. Multi-Point Optimizations. The constrained multi-point optimization has 20 + m design vari- 
ables: the same 20 y-coordinates which describe the geometry and m angles of attack. The minimum angle 
of attack for which the lift constraint is satisfied, decreases with increasing free flow Mach number due to 
compressibility effects [21]. Consequently, each design condition adds one additional angle-of- attack design 
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Fig. 6.4. Drag profile obtained using 4-point optimization strategy (results obtained using Grid 2, Mi = 0.72,0.74,0.76 
and 0.78, w, = 0.25, C ; * = 0.175, and CONMIN optimizer) 

variable. 

The two-point optimization results shown in Figure 6.3 illustrate the shortcomings reported by Drela [9]. 
Optimization at selected Mach numbers results in clearly distinguishable drag troughs at each of the design 
Mach numbers (see also Figures 6.4, 4.2 and 4.3). Drela. leaves it up to the designer to determine which Mach 
numbers to include in the objective in Eq. (4.2) and what weights to choose. Three “reasonable” selections 
are compared with each other in Figure 6.3 and it is clear that, at least in this case, the selection of the 
design conditions has an important effect on the final results. In particular the selection of the endpoints 
of the Mach range can lead to troubling results. We observed this in both 2-point (Figure 6.3) and 4-point 
(Figure 6.4) optimization. 

Since the drag is higher in the upper part of the Mach range, one may want to include more design 
Mach numbers from the upper part than from the lower part in the multi-point objective Eq.(4.2). This 
procedure does not require any additional function calls and consequently keeps the computational cost 
under control and should result in lower drag at the upper end of the Mach range. Figure 6.5 shows the 
result of such a multi-point optimization. The selected Mach numbers are indicated on the chart and the 
weights Wi are obtained from the numerical integration using 4 fixed integration points. It can be concluded 
that the maximum drag is indeed reduced, but a substantial penalty is to be paid near the lower end of the 
Mach range. As a matter of fact the expected value of the drag has actually increased from 0.0044 to 0.0055 
(assuming a uniform distribution for the Mach numbers between 0.7 and 0.8) compared with the result for 
evenly spaced Mach numbers. Also, the distinct multiple drag troughs at each of the sample Mach numbers 
persist, which indicates “localized optimization” as explained earlier. 

All these results indicate that a multi-point optimization does indeed achieve a better overall drag 
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Fig. 6.5. Drag profile obtained using different 4-point optimization strategies (results obtained using Grid 2, C( = 0.6, 
and SLP/ADS optimizer) 

reduction for both methods. This is in line with the findings of other researchers [9]. However, a drag-trough 
is formed at or near each of the discrete design points. The drag increases rapidly away from the design 
points. This is very clear near the high end of the Mach range. 

6.4. Expected Value Optimization. Figure 6.6 shows that the expected value optimization strategy 
results in a much smoother drag profile over the entire Mach range. The integration of the drag over the Mach 
range is performed using spline-based inter /extrapolation. The resulting airfoil geometry is a lot smoother 
as well. It can be concluded that this scheme results in a superior design for an identical computational 
effort. 

For the expected value scheme the integration points in Eq.(5.1) are altered slightly for each optimization 
iteration. Consequently, we can say that the 4-point optimization minimizes a weighted sum of 4 fixed design 
conditions, whereas the robust scheme minimizes the weighted sum of any 4 design conditions. This avoids 
“localized optimization”. In this study a random perturbation scheme is used to change the integration 
points, but other, more adaptive strategies are currently being researched. The number of integration points 
must be sufficiently high so that the integration error caused by the change of integration points between 
optimization steps is smaller than the decrease of drag in that particular optimization step (cfr. supra). 

In his study, Drela [9] applies larger weights to the upper part of the Mach range to ensure “that 
the upper part is not compromised excessively by the less important lower part”. Our statistical decision 
approach suggests that this can only be justified when the higher Mach numbers are more likely to occur, 
i.e. when the PDF /m(M) is concentrated in the upper Mach range. The weights follow automatically from 
the integration procedure and are as such not directly linked to the actual values of the drag Cd(M). 
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Fig. 6.6. Drag profile obtained using different optimization strategies (results obtained using Grid 1 , C( = 0.175 and 
PORT optimizer) 

6.5. Impact of PDF. In most practical cases the Mach number will not be uniformly distributed over 
a given range. A key advantage of the explicitly statistical approach is that the relative importance of each 
operating condition is automatically accounted for through the PDF. 

Quite often a cruise Mach number is set and the assumption of a (truncated) “normal” or Gaussian 
distribution around this mean value seems appropriate. The truncated normal PDF is very well approximated 
by a Beta distribution. The Beta distribution is always bounded within an interval [o, b] and has shape 
parameters a\ and a 2 , both greater than zero: 


/ x — a 

(6.1) Beta(x,a,b,a\,a 2 ) = — — 77- 

where B(a 1,0:2) = , ((*1), ((*2)/, (an + a2) is the Beta-function. 

Depending on the value of ai and a2, the Beta. PDF can assume a variety of shapes (see Figure 6.8). 
When ai = a2 the distribution is symmetric and can be used as an approximation for the truncated normal, 
especially for a* > 5. When ai = a2 = 1, the distribution becomes uniform. For ai ^ a2 the distribution is 
skewed towards either left or right. Bath-tub distributions are obtained when both parameters in the PDF 
are less than 1 (but greater than 0). 

Figure 6.9 compares the optimal drag-profile obtained using three different Beta-distributions bounded 
within [0.7, 0.8]. When a PDF Beta( 3,1) is assumed for the Mach number, the higher Mach numbers are 
more likely to occur. Their impact is easily discerned from Figure 6.9. The greater weight of the higher 
Mach numbers translates into a lower drag in the high Mach range. The trade-off is a higher drag at the 


- a) x B(ai,a2 ) 
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Fig. 6.7. Drag profile obtained using different optimization strategies (results obtained using Grid 2, C( = 0.6 and 
SLP/ADS optimizer) 


lower end of the Mach range. For a Beta(l,3)-distribution, a lower drag is obtained in the lower part of the 
Mach range at the expense of much faster drag increase for higher Mach numbers. When a Beta( 5, 5) PDF 
is used, the higher likelihood of Mach numbers near the mean value M = 0.75 results in the lowest drag 
of all three curves near the middle of the Mach range. The corresponding airfoil geometries are shown in 
Figure 6.10. 

This illustrative comparison reveals the importance of an accurate quantification of the PDF of the 
Mach range. In general, careful data analysis is required when the PDF is selected for each of the uncertain 
variables. Practical experience has revealed that especially the tails of the distributions need to be modeled 
very carefully because they tend to have a large impact [17]. Cut-off values should most definitely not be 
chosen arbitrarily! 

6.6. Approximate Second-Order Second-Moment Method (SOSM). In this section we present 
results using the deterministically equivalent problem formulation Eq.(5.4), which is based on a second-order 
analytic approximation of Eq.(5.1). 

Eq.(5.4) indicates that first-order sensitivities of the drag C& with respect to the uncertain variable M 
do not affect the expected value of the design. The second-order information represents the curvature of the 
Cd{M ) Mach sweep, which is indicative of “localized optimization”. In the SOSM formulation, the weighting 
between the drag and the curvature is determined by physics: it is given by the variance of the Mach number. 
Figure 6.11 shows a considerable improvement in the robustness of the design. 

The drag reduction is not quite as good as obtained using explicit numerical integration (Figure 6.7), 
but the computational effort is a lot smaller. The method is particularly useful if higher order derivatives 
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Probability Density 


o Uniform 



Fig. 6.8. Sample Beta distributions, and comparison with truncated Normal 



Mach number 


Fig. 6.9. Impact of assumed PDF for M on optimal drag profile (grid 2, Cf = 0.6 SLP/ADS optimizer) 
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Solid line: Beta(1 ,3) 
Dashed line: Beta(3,1) 
Dotted line: Beta(5,5) 



Fig. 6.10. Airfoil geometries for three Mach distributions (grid 2, Cf = 0.6 SLP/ADS optimizer) 



Fig. 6.11. Comparison of SOSM result with single-point optimization (grid 2, Cf =0.6 SLP/ADS optimizer) 


are available (and numerically reliable) and several uncertain variables are present in the problem. Table 
6.6 compares the relative computational cost. SOSM scales linearly with the number of random variables, 
whereas full numerical integration may rapidly become prohibitively expensive. Because the method is based 
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Table 6.2 

Number of function/ derivative evaluations required per optimization step (Note: SOSM requires less if analytic derivatives 
are available) 


Optimization Method 

1 Random Variable 

3 Random Variables 

Single Point 

1 

1 

SOSM 

3 

7 

Expected Value 

4 

64 

(using 4-point integration) 




on a second-order approximation of the objective function, SOSM will give the best results if the variance 
of the random variables is relatively small. 

7. Summary and Conclusions. The robustness of an optimal solution can be achieved by incorpo- 
rating the variability of the operating conditions directly into the optimization problem formulation. The 
practical application shows that a statistical approach leads to smooth airfoil geometries and drag profiles. 

The suggested expected value optimization method is computationally similar to existing multi-point 
optimization, which is widely accepted in industry. This increases the likelihood of acceptance by both 
designers and theorists alike. The new formulation avoids the arbitrary selection of design conditions and 
weighting factors; they automatically follow from the procedure. 

A second-order second-moment approximate integration provides additional insight into the problem. 
SOSM scales linearly with the number of design variables and may be the only feasible alternative when a 
lot of uncertainties are considered. This will be the focus of future research. 

The relative likelihood of each operating condition is directly taken into account. The importance of a 
careful selection of the PDF is illustrated by means of an example. A randomized integration scheme ensures 
that the optimizer cannot exploit approximation errors due to discretization. Further research into different 
sampling schemes is required. 

It can be concluded that airfoil shape optimization on the basis of the Euler equations leads to some 
interesting candidate designs. However, viscous effects need to be included to achieve more realistic pressure 
distributions. 
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